TK!2
=v
=u
=l
#1
:n
a0
:v
'a01
'a02
'a03
'a04
'a05
'a06
'a07
'a08
'a09
'a010
E
#2
:n
a01
:v
1
1
1
1
E
#3
:n
a02
:v
1
2
3
4
E
#4
:n
a03
:v
1
3
6
10
E
#5
:n
a04
:v
1
4
10
20
E
#6
:n
a05
:v
E
#7
:n
a06
:v
E
#8
:n
a07
:v
E
#9
:n
a08
:v
E
#10
:n
a09
:v
E
#11
:n
a010
:v
E
#12
:n
a
:v
'a1
'a2
'a3
'a4
'a5
'a6
'a7
'a8
'a9
'a10
E
#13
:n
a1
:v
-.3086863202265788
.787275376005946
-.5303657197721289
.06018672054749688
E
#14
:n
a2
:v
.7230903161940962
-.1632336508926851
-.6403317308849631
.2011726727575942
E
#15
:n
a3
:v
-.594550779583978
-.5321066916667604
-.3918321312058809
.458082328991412
E
#16
:n
a4
:v
.1684117597656997
.2653577325694441
.3938972691795621
.8637521023251448
E
#17
:n
a5
:v
E
#18
:n
a6
:v
E
#19
:n
a7
:v
E
#20
:n
a8
:v
E
#21
:n
a9
:v
E
#22
:n
a10
:v
E
#23
:n
d
:v
.03801601522914029
.4538345500256651
2.203446167647321
26.30470326709787
E
#24
:n
s
:v
1.918509721497501E-27
2.602085213965211E-18
0
0
E
=f
#1
:n
Mcopy
:c
Copy matrix (i.e. subordinate lists)
:t
Procedure
:i
M1,M2
#2
:n
Householder2
:c
Tridiagonalization, eigenvector option
:t
Procedure
:i
a,d,s
#3
:n
QL2
:c
QL algorithm (implicit shifts, eigenvector option)
:t
Procedure
:i
a,d,s
#1
:n
Mcopy
:z
" The procedure copies subordinate lists of the masterlist M1 into 
" subordinate lists of the master list M2.

for i=1 to length(M1)
   call listcopy(M1[i],M2[i])   " error indicates no listname in ith
                                " element of either M1 or M2
next i
/E
E
#2
:n
Householder2
:z
" Notation:  a    on input: original matrix
"                 on output: orthogonal matrix Q (see below)
"            d,s  lists with diagonal and subdiagonal elements of the
"                 transformed matrix

" Description: The procedure reduces the original symmetric matrix to a
" tridiagonal form and generates an orthogonal matrix  Q  such that
" B = PùAùQ  where  A  is the original matrix,  B is the tridiagonal matrix,
" and  P  is the transpose of  Q .  The outputs  a,d,s  are used by the
" procedure QL2 for evaluation of eigenvalues and eigenvectors.

call blank(d)
call blank(s)
n:= length(a[1])
for i=n to 2 step -1
   ii:= i - 1
   h:= 0
   scale:= 0
   if ii<= 1 goto A
     for k=1 to ii
       scale:= scale + abs(a[i][k])
     next k
     if scale=0 goto A
       for k=1 to ii
          a[i][k]:=a[i][k]/scale
          h:= h + a[i][k]^2
       next k
       f:= a[i][ii]
       g:= sqrt(h)
       if f=>0 then g:= -g
       s[i]:= scale*g
       h:= h - f*g
       a[i][ii]:= f - g
       f:= 0
       for j=1 to ii
          a[j][i]:= a[i][j]/h
          g:= 0
          for k=1 to ii
             if k>j then g:= g+a[k][j]*a[i][k] else g:= g+a[j][k]*a[i][k]
          next k
          s[j]:= g/h
          f:= f + s[j]*a[i][j]
       next j
       e:= .5*f/h
       for j=1 to ii
          f:= a[i][j]
          g:= s[j] - e*f
          s[j]:= g
          for k=1 to j
             a[j][k]:= a[j][k] - f*s[k] - g*a[i][k]
          next k
       next j
       goto B
A:   s[i]:= a[i][ii]
B: d[i]:= h
next i
d[1]:= 0
s[1]:= 0
for i=1 to n
   ii:= i - 1
   if d[i]=0 goto C
     for j=1 to ii
        g:= 0
        for k=1 to ii
           g:= g + a[i][k]*a[k][j]
        next k
        for k=1 to ii
           a[k][j]:= a[k][j] - g*a[k][i]
        next k
     next j
C: d[i]:= a[i][i]
   a[i][i]:= 1
   if ii<1 goto D
     for j=1 to ii
        a[i][j]:= 0
        a[j][i]:= 0
     next j
D: next i
/E
E
#3
:n
QL2
:z
" Notation:  a    on input:  orthogonal matrix Q (see function Householder2)
"                 on output: a matrix with columns representing eigenvectors
"            d,s  on input:  lists with diagonal and subdiagonal elements
"                            of a tridiagonal matrix
"                 on output: list  d  contains eigenvalues of the original
"                            matrix

" Description: The procedure performs a sequence of orthogonal transfor-
" mations yielding the eigenvalues and eigenvectors of a given symmetric
" matrix.  Loss of accuracy is prevented by implicit shifts.

n:= length(d)
call listcopy(s,s,2)
s[n]:= 0
for k=1 to n
   iter:= 0
A: for m=k to n-1
      D:= abs(d[m]) + abs(d[m+1])
      if (abs(s[m])+D)=D goto B
   next m
   m:= n
B: if m=k goto E
     if iter=30 then error:= too_many_iterations
     iter:= iter + 1
     g:= (d[k+1]-d[k])/(2*s[k])
     r:= sqrt(g^2+1)
     if g<0 then r:= -r
     g:= d[m] - d[k] + s[k]/(g+r)
     w:= 1
     c:= 1
     p:= 0
     for i=m-1 to k step -1
        f:= w*s[i]
        b:= c*s[i]
        if abs(f)<abs(g) goto C
          c:= g/f
          r:= sqrt(c^2+1)
          s[i+1]:= f*r
          w:= 1/r
          c:= c*w
          goto D
C:        w:= f/g
          r:= sqrt(w^2+1)
          s[i+1]:= g*r
          c:= 1/r
          w:= w*c
D:      g:= d[i+1] - p
        r:= (d[i]-g)*w + 2*c*b
        p:= w*r
        d[i+1]:= g + p
        g:= c*r - b
        for j=1 to n
           f:= a[j][i+1]
           a[j][i+1]:= w*a[j][i] + c*f
           a[j][i]:= c*a[j][i] - w*f
        next j
     next i
     d[k]:= d[k] - p
     s[k]:= g
     s[m]:= 0
     goto A
E: next k
/E
E
=r
#1
:r
" ***  Eigenvalues and Eigenvectors of Symmetric Matrices  ***
:s
C
#3
:r
call Mcopy('a0,'a)
#4
:r
call Householder2('a,'d,'s)
#5
:r
call QL2('a,'d,'s)
#7
:r
" Notes:
:s
C
#9
:r
" This model illustrates evaluation of eigenvalues and eigenvectors of
:s
C
#10
:r
" symmetric matrices following algorithms implemented in the EISPACK package
:s
C
#11
:r
" and referred to in "Numerical Recipes" by W.H.Press et al. (Cambridge
:s
C
#12
:r
" University Press, 1986).
:s
C
#14
:r
" The original matrix is represented by a master list  a0  and subordinate
:s
C
#15
:r
" lists a01, a02, ..., which can be set up in the interactive table
:s
C
#16
:r
" matrix_a0 .  Upon completing the call to the function  Householder2  the
:s
C
#17
:r
" lists  d  and  s  hold the diagonal and subdiagonal coefficients of the
:s
C
#18
:r
" tridiagonalized matrix (they may be seen in the table  ds ).  The matrix 
:s
C
#19
:r
" a  with the subordinate lists  a1, a2, ...  in table  matrix_a  contains
:s
C
#20
:r
" the matrix  Q  used later by the function  QL  for evaluation of the 
:s
C
#21
:r
" eigenvectors.
:s
C
#23
:r
" After a call to  QL2  is completed, list  d  contains the eigenvalues, and
:s
C
#24
:r
" the  kth  column of the matrix  a  represents a normalized eigenvector 
:s
C
#25
:r
" corresponding to the  kth  eigenvalue.
:s
C
#27
:r
" The functions  Householder2  and  QL2  may be run in a row or separately
:s
C
#28
:r
" (cancelling first the 3rd rule, and then the first two).  The generation
:s
C
#29
:r
" and transformation of the matrices and lists may be watched in the inter-
:s
C
#30
:r
" active tables.
:s
C
=p
=t
#1
:n
matrix_a0
:s
Screen
:v
Vertical
:f
1
:r
""
:c
""
:t
"Original symmetric matrix"
:u
#1
:l
a01
:w
6
:h
""
#2
:l
a02
:w
6
:h
""
#3
:l
a03
:w
6
:h
""
#4
:l
a04
:w
6
:h
""
#5
:l
a05
:w
6
:h
""
#6
:l
a06
:w
6
:h
""
#7
:l
a07
:w
6
:h
""
#8
:l
a08
:w
6
:h
""
#9
:l
a09
:w
6
:h
""
#10
:l
a010
:w
6
:h
""
E
#2
:n
matrix_a
:s
Screen
:v
Horizontal
:f
1
:r
""
:c
""
:t
"Original matrix turning into elements of Q matrix and eigenvectors"
:u
#1
:l
a1
:w
6
:h
""
#2
:l
a2
:w
6
:h
""
#3
:l
a3
:w
6
:h
""
#4
:l
a4
:w
6
:h
""
#5
:l
a5
:w
6
:h
""
#6
:l
a6
:w
6
:h
""
#7
:l
a7
:w
6
:h
""
#8
:l
a8
:w
6
:h
""
#9
:w
10
:h
""
E
#3
:n
td
:s
Screen
:v
Vertical
:f
1
:r
""
:c
""
:t
"Diagonal and subdiagonal elements turning into eigenvalues and 0's"
:u
#2
:l
d
:w
10
:h
""
#4
:l
s
:w
10
:h
""
#6
:w
10
:h
""
E
=n
=g
#1
:c
.000001
:i
10
:v
Yes
:l
60
:f
66
:m
0
:w
80
:b
Yes
:a
Yes
:r
Yes
:h
No
:t
1
:n
Yes
:d
LPT1
%Tr,f,0,1,0,1,0
